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^"*^ . Abstract. Discrete- velocity approximations represent a popular way for computing the 

Boltzmann collision operator. The direct numerical evaluation of such methods involve 
a prohibitive cost, typically 0(N 2d+1 ) where d is the dimension of the velocity space. In 
this paper, following the ideas introduced in [271128] . we derive fast summation techniques 
for the evaluation of discrete- velocity schemes which permits to reduce the computational 
cost from 0(N 2d+1 ) to 0{N d N d log 2 N), N < N, with almost no loss of accuracy. 
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1. Introduction 

Among deterministic methods to approximate the Boltzmann collision integral, one of 
the most popular is represented by discrete velocity models (DVM). These methods [7, 
ESI El EH EH1 E3 123 E] are based on a regular grid in the velocity field and construct a 

l 
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discrete collision mechanics on the points of the quadrature rule in order to preserve the 
main physical properties. 

As compared to Monte-Carlo methods, these methods have certain number of assets: 
accuracy, absence of statistical fluctuations, and the fact that the distribution function 
is explicitly represented in the velocity space. However their computational cost is more 
than quadratic and they cannot compete with the linear cost of a Monte Carlo approach. 
Indeed the "naive" cost of a product quadrature formula for the (d—l) + d fold Boltzmann 
collision integral in dimension d is 0(M d ~ 1 N d ), where M is related to the angle and TV" 
to the velocity discretizations. More concretely Buet presented in [7] a DVM algorithm 
widely used since then in 0(N 2d+1+e ) for all e > (and a constant depending on e); 
Michel and Schneider algorithm in [26J is 0(N 2d+s ) where 5 depends on d and is close 
to 1; finally the method of Panferov and Heinz |3U] is 0(N 2d+l ). For this reason several 
acceleration techniques for DVM have been proposed in the past literature. We do not 
seek to review them here, and refer the reader to [22j EU 155], 157]. 

More recently a new class of methods based on the use of spectral techniques in the 
velocity space has attracted the attention of the scientific community. The method first 
developed for the Boltzmann equation in [3T] is based on a Fourier- Galer kin approxima- 
tion of the integral collision operator. As shown in [321 [33] the method permits to obtain 
spectrally accurate solution at a reduced computational cost of 0(N 2d ). A proof of sta- 
bility and convergence for this method has been given in [16]. Finally the method has 
been extended to the case of the quantum Boltzmann collision operator [151 120] . Other 
methods based on spectral techniques have been developed in [31 [T5]. 

One of the major differences between DVM and spectral methods is that in the latter 
the interaction kernel of the Boltzmann collision integral is not modified in order to obtain 
a conservative equation on a bounded domain. This aspect has a profound influence on 
the resulting structure of the algorithm since most of the symmetries which are present in 
the original operator are preserved. Using this fact, in [27J [2H], the authors developed a 
numerical technique based on the Fast Fourier Transform (FFT) that permits to reduce 
the cost of spectral method from 0(N 2d ) to 0(M d ~ l N d log 2 N) where M is the number 
of angle discretizations. These ideas have been successfully used in [T7] to compute space 
non homogeneous solutions of the Boltzmann equation. 

In this paper we will consider general discrete velocity approximation of the Boltzmann 
equation without any modification to the original collision kernel and show how the FFT 
techniques developed in j27J [2H] can be adapted to this case to obtain acceleration algo- 
rithms. In this way, for a particular class interactions using a Carleman-like representation 
of the collision operator we are able to derive discrete velocity approximations that can 
be evaluated through fast algorithms at a cost of 0(N d N d log 2 N), N <C N. The class of 
interactions includes Maxwellian molecules in dimension two and hard spheres molecules 
in dimension three. 

Let us emphasize here that a detailed analysis of the computational complexity in DVM 
is non trivial since imposing conservations on the points of the quadrature rule originates 
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a summation formula that requires the exact enumeration of the set of involved orthogonal 
directions in 

The rest of the paper is organized in the following way. In the next Section we introduce 
briefly the Boltzmann equation and give a Carleman-like representation of the collision 
operator which is used as a starting point for the development of our methods. In Section [3] 
a fast DVM method is introduced together with a detailed analysis of its computational 
complexity. In Section U we present some numerical results obtained with the fast and 
the classical DVM methods. 



2. Preliminaries 

2.1. The Boltzmann equation. The Boltzmann equation describes the behavior of a 
dilute gas of particles when the only interactions taken into account are binary elastic 
collisions. It reads for x, v € M rf (d> 2) 

^ + v-V x f = Q(f,f) 

where f(t,x,v) is the time-dependent particle distribution function in the phase space. 
The Boltzmann collision operator Q is a quadratic operator local in (t,x). The time and 
position acts only as parameters in Q and therefore will be omitted in its description 

(2.1) Q(f,f)(v)= [ f B{cos9,\v-v*\) [f'j'-fj] da dv*. 

In (12. ip we used the shorthand / = /(«),/* = /(f*)> / = f( v ')i f* = f( v '*)- The velocities 
of the colliding pairs (v,v*) and (v',vl) are related by 

, v + v* \v — vJ . v + v* \v — vJ 

v = 1 a, v ± = a. 

2 2 ' * 2 2 

The collision kernel B is a non-negative function which by physical arguments of invariance 
only depends on \v — v+\ and cos 9 = g ■ a (where g = (v — — 

Boltzmann's collision operator has the fundamental properties of conserving mass, mo- 
mentum and energy 

Q(f,f)<Kv)dv = Q, <t>(v) = l,v,\v\ 2 
and satisfies the well-known Boltzmann's i?-theorem 

~ [ flogfdv = -[ Q(/,/)log(/)cfo>0. 

at Jv£R d Jv£K. d 

The functional — J f log / is the entropy of the solution. Boltzmann i?-theorem implies 
that any equilibrium distribution function, i.e. any function which is a maximum of the 
entropy, has the form of a locally Maxwellian distribution 

M( P ,u,T)(v) = -^JL-p exp f-^M , 
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where p, u, T are the density, mean velocity and temperature of the gas 

p= f f(v)dv, u=~ [ vf(v)dv, T = — ( \u — v\ 2 f(v)dv. 

JvGR d P Jv€R d dp J V £R d 

For further details on the physical background and derivation of the Boltzmann equation 
we refer to and [3S]- 

2.2. Carleman-like representation in bounded domains. In this short paragraph 
we shall approximate the collision operator on a bounded domain starting from a repre- 
sentation which somehow conserves more symmetries of the collision operator when one 
truncates it in a bounded domain. This representation was used in [H El [211 12H] and is 
close to the classical Carleman representation (cf. |10j). 
The starting point of this representation is the identity 

(2.2) - / F(\u\a - u) da = , ,, „ / 5(2 x ■ u + \x\ 2 ) Fix) dx. 

2 J$d-1 \u\ a Z Ju d 

It can be verified easily by completing the square in the delta Dirac function, taking the 
spherical coordinate x = ra and performing the change of variable r 2 = s. Then, setting 
u = v — i>* and r = |it|, we have the following Lemma. 

Lemma 2.1 (Cf. [28], subsection 2.1). Introducing the change of variables 

1 

x = - ra, y = — v — x, 
the collision operator 112.1]) can be rewritten in the form 

Q(f, /)(«) = / / B(x, y) 8(x ■ y) [f(v + y) f(v + x )-f( v + x + y ) /(„)] dx dy, 

Jx£R d Jy£R d 

where 

(2.3) B(x,y) = B(\x\, \y\) = 2^ B ( ]*\ J\x\* + \yA (\x\ 2 + |y| 2 )"^. 

Now let us consider the bounded domain T>t = [—T,T] d (0 < T < +oo). First one 
can remove the collisions connecting with some points out of the box. This is the natural 
preliminary stage for deriving conservative schemes based on the discretization of the 
velocity. In this case there is no need for a truncation on the modulus of x and y since we 
impose them to stay in the box. It yields 

Q tr (fJ)(v)= I i A B(x,y)S(x-y) 

J J{x,yeR d | v+x, v+y, v+x+y G X> T j 

[f(v + y) f(v + x) - f(v + x + y) f(v)] dx dy 
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defined for v £ T>t- One can easily check that the following weak form is satisfied by this 
operator 

(2.4) I Q tT {fJ)<f{v)dv = - f ( f B(x,y)S(x-y) 

J Q J J J \v,x,y £R d | v, v+x,v+y, v+x+y £T> T j 

f(v + x + y) f{y) [(p(v + y) + ip(v + x) — ip(v + x + y) — <p(v)] dv dx dy 

and this implies conservation of mass, momentum and energy as well as the //-theorem on 
the entropy. The problem of this truncation on a bounded domain is the fact that we have 
changed the collision kernel itself by adding some artificial dependence on v, v*, v' , v'^. In 
this way convolution-like properties are broken. 

A different approach consists in truncating the integration in x and y by setting them 
to vary in Br, the ball of center and radius R. For a compactly supported function / 
with support Bs, we take R = S in order to obtain all possible collisions. Since we aim at 
using the FFT algorithm to evaluate the resulting quadrature approximation, and hence 
we will make use of periodic distribution functions, we must take into account the aliasing 
effect due to periods superposition in the Fourier space. As for the spectral method a 
geometrical argument (see [32] for further details) shows that using the periodicity of the 
function it is enough to take T > (3 + y/2)S/2 to prevent intersections of the regions where 
/ is different from zero. 

The operator now reads 

(2.5) Q R (f,f)(v)= [ I B{x,y)5{x-y) 

[f(v + y)f(v + x) - f(v + x + y)f(v)] dx dy 

for v S T>t- The interest of this representation is to preserve the real collision kernel and 
its properties. It is easy to check that, except for the aliasing effect, the operator preserves 
all the original conservation properties, see the weak form in equation (j2.6|) . 

In order to understand the possible effect of periods superposition we can rely on the 
following weak form valid for any function ip periodic on Vt 

(2.6) / Q R (f,f)<p(v)dv = - -f { ( B(x,y)6(x-y) 
JV T 4 Jv£V T Jx&Br JyeBR 

f(v + x + y)f{v) [ip(v + y) + ip(v + x) — ip(v + x + y) - f(v)\ dv dx dy. 

About the conservation properties one can show that 

(1) The only invariant ip is 1: it is the only periodic function on T>t such that 

<p(v + y) + (f(v + x) - ip(v + x + y) - <p(v) = 

for any v £ T>t and X-Ly £ Br (see jll] for instance). It means that the mass is 
locally conserved but not necessarily the momentum and energy. 

(2) When / is even there is global conservation of momentum, which is in this case. 
Indeed Q R preserves the parity property of the solution, which can be checked 
using the change of variable x — > — x, y — > — y. 
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(3) The collision operator satisfies formally the i?-theorem 

/ Q R (f,f)log(f)dv<0. 

(4) If / has compact support included in 3$ with T > (3 + y/2)S/2 (no-aliasing con- 
dition, see [32] for a detailed discussion) and R = S, then no unphysical collisions 
occur and thus mass, momentum and energy are preserved. Obviously this com- 
pactness is not preserved with time since the collision operator spreads the support 
of / by a factor y/2. 

2.3. Application to discrete-velocity models. The representation Q R of this section 
can also be used to derive discrete velocity models (DVM). Any DVM can be written in 
the general form 

(2-7) A(/,/) = E 'tjA/} !,f L . 

j,k,iez d 

where Di denotes the discrete Boltzmann collision operator and the integer indexes refer 
to the points in the computational grid. 

In order to keep conservations the coefficients T^j are defined by 

(2.8) rf j = 1(< + j - k - I) l(|i| 2 + |i| 2 - |Af - |/| 2 ) B(\k -i\,\l- j\) ti# 

where 1 denotes the function on Z defined by l(z) = 1 if z = and elsewhere, and 
w i 'j > are the weights of the quadrature formula, which characterize the different DVM. 
The function B > is the discrete collision kernel. One can check on this formulation 
that the scheme satisfies the usual conservation laws and entropy inequality (see [Ml [S] 
and the references therein). More details on the DVM schemes can also be found in [B]. 

Thanks to equations (|2.7|) and (|2.8|) . we can write at the discrete level the same repre- 
sentation as in the continuous case 

A (/>/)= X/ ^k,l[fi+kfi+l ~ fifi+k+l] 

with 

f k,i = 2^ B (^=pr V^ + I'l 2 ) + I'l 2 )-^ • 

This is coherent with the DVM obtained by quadrature starting from the Carleman rep- 
resentation in |3U| . 

Now again when one is interested to compute the DVM in a bounded domain there are 
two possibilities. First as in the case of Q^ T one can force the discrete velocities to stay in 
a box, which yields for i G [—A, N~l d (again using the one index notation for d-dimensional 
sums) 

D Y (/; /) = ^2 Tfc,z [fi+kfi+i - fifi+k+i] ■ 
k,i 

—N< i+k, i+l, i+k+l<N 
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This new discrete operator is completely conservative but the collision kernel is not invari- 
ant anymore according to i, which breaks the convolution properties and then prevents 
the derivation of a fast algorithm. 

The other possibility is to periodize the function / over the box and truncate the sum 
in k and I. It yields for a given truncation parameter A'gN* 

( 2 -9) Df (/,/)= Y ^k,i[fi+kfi+i - fifi+k+i], 

-N<k,l<N 

for any ie{-N,N} d . 

It is easy to see that D N satisfies exactly a discrete weak form and conservation pro- 
perties similar to Q R . Let us briefly state and sketch the proof of the conservation and 
stability properties of the scheme. 

k I 

Proposition 2.2. Assume that the quadrature weight w i 'j > are positive. Consider 
some truncation numbers N < N G N* and some non-negative initial data fi(0) > 0, 
i G [— N, N~l d . Then the discrete evolution equation 

dtfi = D?(f,f)= Y Fk,l[fi+kfi+i-fifi+k+i], iel-N,Nf, 

-N<k,l<N 

is globally well-posed in M.^~ N > N ^ d . Moreover the coefficients fi(t) are non-negative for all 
time, and 

vt>o, Y f(t)= Y W)- 

iel-N,N] d iel-N,N] d 

Remark 2.3. The DVM scheme we consider therefore preserves non-negativity, but let 
us also emphasize that it preserves momentum and energy up to aliasing issues. This 
is different from spectral methods where the truncation of Fourier modes introduces an 
additional error in the conservation laws. Concerning the spectral method, stability and 
convergence have been proved recently in [T7j to hold in L 1 but only asymptotically, i.e. 
for N big enough related to the initial data. 

Proof of Proposition \2.2i We have the following L 1 -like estimate 
d 



E W)\= E 

ie[-N,Nj d iel-N,Nj d 



E ^k,l[fi+kfi+l ~ fifi+k+l] 



-N<k,l<N 
2 



(2-10) < C Y \fi 

\ie{-N,NY 

The use of a Gronwall argument then gives the local well-posedness of the scheme in 
jp-NjiV] _ Moreover, given a local solution fi(t), for t £ [0,T] and T > 0, it is clear by 
construction that the conservation of mass holds. 

The proof of preservation of non-negativity for this solution is essentially contained in 
the pioneering work of Carleman [ID]. We will sketch its proof in the following. Let us 
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rewrite the system of ordinary differential equations satisfied by fi for a fixed i € [— N, Nj d 

as 

(2-11) ^7 fi + fi E fi+k+l= E ^k,l fi+k fi+l- 

" -N<k,l<N -N<k,l<N 

Let us assume by contradiction that we have 

|7;(t)>0, Vt€[0,T[, VjG l-N,Nf, 

\fi(T) = 0. 
Then, we have necessarily 

fl(T)<0, 

and thus, according to (|2,11|) . 

E r k ,lfi+k(T)fi+i(T)<0. 

-N<k,l<N 

By continuity in time of fj, it comes that 

/ i (T) = 0, VjG [-iV,iV] d 

As these conditions implies using (|2.1ip that = for all t € [0,T], we have a 

contradiction with the non-negativity of the initial condition. 

Finally, the conservations of mass and non-negativity implies the preservation of L 1 
norm, and we can iterate the argument giving the local well-posedness (still using inequal- 
ity (|2.10p ) to obtain the global well-posedness of the scheme. 

□ 

Finally one can derive the following consistency result from [Si)\ Theorem 3] in the case 
of hard spheres collision kernel with d = 3 

Theorem 2.4. Assume that f,g£ C fc (IR 3 ) (k>l) with compact support Bs- The uniform 
grid of step h is constructed on the box T>t with the no-aliasing condition T > (3+v / 2)>S'/2. 
Then for N = [S/h] (where [•] denotes the floor function) and h > sufficiently small, 

\Q(g,f)-D*(g,f) 

where D N is the DVM operator defined in (|2.9|) (for the precise quadrature weights derived 
in [30]J on the grid above-mentioned, and fi = f(ih). Here r = k/(k + 3) and the constant 
C is independent on h. 

Remark 2.5. As can be seen from Theorem 12. 4| the periodized DVM presented in 
this subsection is expected to have a quite poor accuracy. On the contrary the spectral 
method [51] . even in the fast version of [28J, has been proven to be spectrally accurate, 
i.e. of infinite order for smooth solutions. Nevertheless this periodized DVM has some 
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interesting features compared to the spectral method: preservation of sign, stability, and 
preservation of the conservation laws up to aliasing issues. 

3. Fast DVM's algorithms 

The fast algorithms developed for the spectral method in [28] can be in fact extended 
to the periodized DVM method. The method that originates was triggered by the reading 
of the direct FFT approach proposed in [H El H] . 

3.1. Principle of the method: a pseudo-spectral viewpoint. We start from the 
periodized DVM in [— N, N^ d with representation (|2.9|) and as in the continuous case we 
set, for k, I £ -N <k,l<N, 



B(\k\,\l\) = 3 d - 1 B [ ^ +w i y/\W + \i\ 2 ) (l*l 2 + \l\ 2 r^- 

With this notation 

f M = l(k-l)B{\k\,\l\)w k;l , 

and thus the DVM becomes 

dtfi= E • I) B(\k\, \l\)w kj i [fi+kfi+i - fifi+k+l]- 

-N<k,l<N 

Now we transform this set of ordinary differential equations into a new one using the 

involution transformation of the discrete Fourier transform on the vector (fi)—N<i<N- 

This involution reads for I £ [— N, Nj d 

I N N 

ft = 9AT + 1 E /i e -/W> /»= E flZltt) 

^ i=-N I=-N 

where e^(/c) denotes e 2N + 1 , and thus the set of differential equations becomes 

N ( I N 
d *fi = E ^TTT E e K+L-i( 

K,L=-N \ T i=-N 

l(k ■ I) B(\k\, \l\) w Kl (e K (k)e L (l) - e L (k + /)) ,/ /v f L 

-N<k,l<N 

for —N<I<N. We have the following identity 

1 N 

^— £ e K+L ^(i) = 1(K + L - I) 

i=—N 

and so the set of equations is 

JV 

(3-1) d t h= E P(K,L)f K f L 

K,L=-N 
K+L=I 
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with 

P(K,L)= ]T l(k-l)B(\k\,\l\)w k j[e K (k)e L (l)-e L (k + l)} = (3(K, L) — (3(L, L) 

-N<k,l<N 

where 

(3.2) P(K,L)= Yl Hk-l)B(\k\,\l\)w k je K (k)e L (l). 

-N<k,l<N 

Let us first remark that this new formulation allows to reduce the usual cost of computation 
of a DVM exactly to 0(N 2d ) (as with the usual spectral method) instead of 0(N 2d+5 ) 
for 5 ~ 1 [U Ell ED]. Note however that the (2N + l) d x (2N + l) d matrix of coefficients 
(f3(K, L))k,l has to be computed and stored first, thus the storage requirements are larger 
with respect to usual DVM. Nevertheless symmetries in the matrix can substantially reduce 
this cost. 

Now the aim is to give an expansion of (3(K, L) of the form 

M 

(3 K , L ^J2a p (K)a' p (L), 
P =l 

for a parameter M £ N* to be defined later. Indeed, this formulation will allow us to 
write (|3.1|) as a sum of M discrete convolutions and then this algorithm can be computed 
in 0(M N d log 2 (-/V)) operations by using standard FFT techniques [El [9], as in the fast 
spectral method. 

3.2. Expansion of the discrete kernel modes. We make a decoupling assumption on 
the collision kernel as in the spectral case [25] 

(3.3) B(\k\,\l\)w k)l = a(k)b(l). 

Note that the DVM constructed by quadrature in dimension 3 for hard spheres in [3D] on 
the cartesian velocity grid h Z 3 (for h > 0) satisfies this decoupling assumption with a{k) = 
h 5 \k\/gcd(k\, k2, /C3) and b(l) = 1 (see [SUJ Formula (20-21)]), and gcd(/ci, &2, £3) denotes 
the greater common divisor of the three integers. For Maxwell molecules in dimension 2 
on the grid /1Z 2 , these coefficients are a{k) = h? \k\/gcd{ki, /C2) and b(l) = 1. 

The difference here with the spectral method, which is a continuous numerical method, 
is that we have to enumerate the set of {— N < k,l < N \ k.Ll}. This motivates for a 
detailed study of the number of lines passing through and another point in the grid (this 
is equivalent to the study of this set), in order to compute the complexity of the method 
in term of N. 

To this purpose let us introduce the Farey series and a new parameter < N < N for 
the size of the grid used to compute the number of directions. The usual Farey series is 

f% = {(J>,<l) G [0, Nj 2 \ < p < q < N, q>l, and gcd{p,q) = l} 

where gcd(p, q) denotes again the greater common divisor of the two integers (more details 
can be found in [19J). We gave a schematic representation of the two dimensional Farey 
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Figure 1. Representation of the Farey series and of A 1 ^, the primal 
representant of lines in [-N, N] associated, for N = 7 and N = 3 



series in Figure [TJ It is straightforward to see that the number of lines A 1 ^ passing through 
in the grid \-N, N] 2 is 



A 1 - 

N 



( F 1 - 

\ N 



1 



where the factor 4 allows to take into account the permutations when counting the couples 
(j),q) as well as the ordering (symmetries in Figure 1), minus the line which is repeated 
during the symmetry process. 
Similarly one can define the set 



F% = {(P,Q,r) G [0, iVj 3 | 0<p<q<r<N, r > 1, and gcd{p,q,r) = l} 
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and the number of lines A 2 ^ passing through in the grid [— N, iV] 3 is 



A 



N 



24 



P 

J i 



N 



T 1 - 

N 



2^ 



all possible permutations of the three numbers times 4 and minus the interfaces 2A\ 



N 



accounting for the possible negative values by symmetry, minus 24 LFl for the spurious 
terms when two equal numbers are swapped. The exponents of the Farey series refer to 
the dimension of the space of lines (which is d— 1). Now let us estimate the cardinals of 
JFi and J^. 



Lemma 3.1. The Farey series in dimension d 
asymptotic behavior 



2 and d = 3 satisfy the following 



F 1 - 

N 



N 2 



N 



2<(2) 
N 3 



+ 0(N logiV) 



+ 0(N 2 



3N 2 

— + 0(N logiV), 

7T Z 



12C(3) 

where £(s) = 2~) n >o n ~ s denotes the usual Riemann zeta function. 
Remark 3.2. In dimension d, the formula would be 
Fff 1 = {(pi,P2, ■ ■ ■ ,Pd) € lO,Nf | < pi < P 2 < ■ ■ ■ < p d < N, Pd > 1 

and gcd(pi,p 2 , 

The cardinal of Fff 1 could be computed by induction with the same tools as in the proof: 



■d-l 

N 



C d — + 0(N d - 1 ) 



((d) 



1 



The non-negative constant is given by 

Cd ' r ~ 2 d ~ 2 d\ ' 

the factorial coming from the successive summations of the Riemann series. 

Proof of Lemma \S.1\ The proof of the first equality is extracted from [191 Theorems 330 
& 331 page 268], and given shortly for convenience of the reader. The proof of the second 
inequality is inspired from this first proof. 

Let us introduce (p(n) the Euler function (i.e. the number of integers less than and 
prime to n) and the multiplicative Mobius function /i(n) such that /i(l) = 1, fJ>(n) = 
if n has a squared factor and ^(p\P2 ■ • • Pk) = ( — if all the primes pi,P2, ■ ■ ■ ,Pk are 
different. We have the following connection between these two arithmetical functions (see 
[HI Formula (16.3.1), page 235]): 

Kd) 



cp{n) 



nJ2 

d\n 



d 



dd'=n 
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Now let us compute the cardinal of the Farey series in dimension 2: 
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N 



•'n 



<p(i)+...+ip(N) = e e d '^d) 

m=l dd'=m 
N ([N/d] 

E = 1 + E m(<o E rf/ 



dd'<N 
N 



d=l \d'=l 



1 JV - TV 



2 d =l ' 7 2 



TV 1 \ at2 oo / j\ / oo 



- o(N log iV 



<Z=1 
iV 2 
2C(2) 



2 £ d2 



d 



+ 0{N) + 0[N\ogN 



N+l 



d 2 



where we have used the classical formula l/C( s ) = X^Li n(n)/n s (cf. [THl Theorem 287, 
page 250]). 

Now for the dimension d = 3, we enumerate the set J-^ in the following way: we fix 
r > 1 then 1 < q < r (the case q = is trivial and treated separately), then p < q such 
that gcd(p, gcd(g, r)) = 1 (we use the associativity of the function gcd). This leads us 
to count the number of p in such that gcd(p,5) = 1 for a given <5|g. When <5 > 1, 

writing p = k5 +po with po G [[1, <5 — 1]] , this number is seen to be ip(5) (q/S). When 5=1 
this number is q + 1 (all the values from to q). Thus the formula (p(5) (q/S) is still valid 



if we deal separately with the case p = 0, which has cardinal 
the cardinal of J-^ . We first write 



N 



Now let us compute 



(3.4) 



N 



N r 



r=l q=l 



y?(gcd(g,r)) 
gcd(g, r) 



r = lq=l d \q,d\r ° 
i N t j\ N r 

d=l r=lq=l 
d\r d\q 
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We shall now focus on the last member of the right hand side of this expression. We have 



N r 



N [r/d] 



d 



N /r.i2 



EE« = <*EE* = *E 

-1 d! 

d\r 

[N/d] 



r=l 5=1 
d\r d\q 



r=l d'=l 
d\r 



r=l 
d\r 



+ 



(3.5) 



Finally, we obtain by plugin (|3.5p into ([37 



f E ((<n 2 + d" 

d"=l 

\ (\ (N/df + O ((N/d) 2 ) + O (N/d)) . 



iV 



O \ N ) + 



1 JV /I _ 3 

- £ M (c0 - (iV/d) + O ((N/df) + O (N/d) 



d=l 



d=l 



2 i 

N 3 



1 



d=N+l 



" 12C(3) 
This conclude the proof. 



+ 



O (n 2 ) 



□ 



Now one can deduce the following decomposition of the kernel modes using their defi- 
nition (|3.2|) and the decoupling assumption (|3.3|) on the discrete kernel 



P(K,L) = J2 Hk.l)a(\k\)b(\l\)e K (k)e L (l) 

-N<k,l<N 

* f(K,L)= J2 [ E a(\k\)e K (k)}[ £ b(\l\) e L (l) 



N —N<k<N 



-N<KN 



with equality if N = N. Here 1 denotes the set of primal representants of directions 
of lines in [— N, NJ passing through 0. After indexing this set, which has cardinal , 
one gets 

(3.6) 



P*(K,L)= a P (K)a' p (L) 
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with 

a p (K)= a' p (L)= £ b(\l\)e L (l). 

k£e p Z ; ee J- 
-N<k<N -N<l<N 

After inversion of the discrete Fourier transform, this method yields a decomposition of 
the discrete collision operator 

A d-i 

(3.7) Df ~ D?> N = ]T D?'"*, i G l-N, Nf, 

P =i 

with equality with (|2.9|) if N = N. Each D^' N,p (f, f) is defined by the p-th term of 

the decomposition of the kernel modes (13. 6|) . Each term D N,N ' P of the sum is a discrete 
convolution operator when it is written in Fourier space. Moreover, each a p (resp. a' p ) is 
defined as the discrete Fourier transform of some non- negative coefficients a(|fc|) times the 
characteristic function of k G e p Z (resp. b(\l\) times the characteristic function of I G e p ). 

Hence, we get after inversion of the transform that D N,N ' P is a discrete convolution with 
non-negative coefficients. 

By using the approximate kernel modes j5 N (K, L) , we obtain a new discrete evolution 
equation, which heritates the same nice stability properties as the usual DVM schemes, 
as stated in the following proposition. Its proof is exactly similar to the one of Proposi- 
tion 12.21 when computing by inverse Fourier transform the coefficients associated to 

the approximate kernel modes j3 N (K,L). 

Proposition 3.3. Assume that the quadrature weight w^j > are positive. Consider 
some truncation numbers N < N < N G N* and some non-negative initial data /i(0) > 0, 
% G J— N, N~l d . Then the discrete evolution equation 

(3.8) dtfi = D?>*(f,f), iel-N,N] d , 

is globally well-posed in M.^~ N,N ^ d . Moreover the coefficients fi(t) are non-negative for all 
time, and 

w>o, Yl E /*(°)- 

iel-N,N} d iel-N,N] d 

Remark 3.4. Using the non-negativity of the coefficients together with the conservation 
of mass, momentum and energy, we can prove thanks to standard arguments (see (Hjj that 
the discrete entropy of solutions to the fast DVM method is non-increasing in time. 

3.3. Implementation of the algorithm. The fast DVM method described in the last 
subsection depends on the three parameters iV (the size of the gridbox), R (the truncation 
parameter) and N (the size of the box in the space of lines). The only constraint on these 
parameters is the no-aliasing condition that relates R and the size of the box (and thus R 
and N, thanks to the parameter N). 
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Thus one can see thanks to Lemma 13.11 that even if we take N = N = N, i.e. we take 
all possible directions in the grid [— N, Nj d , we get the computational cost 0(N 2d log 2 N) 
which is better than the usual cost of the DVM, 0(N 2d+1 ) (but slightly worse than the 
cost 0(N 2d ) obtained by solving directly the pseudo-spectral scheme, thanks to a bigger 
storage requirement). 

More generally for a choice of N < N we obtain the cost 0(N d N d log 2 N), which is 
slightly worse than the cost of the fast spectral algorithm (namely 0(M d ~ 1 N d log 2 N) 
where M is the number of discrete angle [28]), but interesting given that the algorithm 
is accurate for small values of N, and more stable. The justification for this is the low 
accuracy of the method (the reduction of the number of direction has a small effect on the 
overall accuracy of the scheme) . 

Finally, as for the fast spectral algorithm, the decomposition (|3.7|) is completely paral- 
lelizable and the computational cost should be reduced (formally) on a parallel machine 
up to 0(N d log 2 N). This method also has the same adaptivity of the fast spectral algo- 
rithm: in a space inhomogeneous setting, the parameter N can be made space dependent, 
according to the fact that some regions in space deserve less accuracy than others, being 
close to equilibrium. 
Remark 3.5. 

(1) Concerning the construction of the set of directions ■A'j^, it can be done with sys- 
tematic algorithms of iterated subdivisions of a simplex, thanks to the properties of 
the Farey series. In dimension d = 2 this construction is quite simple (see |19j ). 
In dimension 3 we refer to [2JJ] . 

(2) Let us remark that in order to get a regular scheme (i.e with no other conservation 
laws than the usual ones) in spite of the reduction of directions, it is enough that 
the schemes contains the directions and tt/2 (see This is satisfied if we 
take the directions contained in J^f -1 , i.e. as soon as N > 1. 

(3) Finally in the practical implementation of the algorithm one has to take advan- 
tage of the symmetry of the decomposition (13. 6h in order to reduce the number of 
terms in the sum: for instance in dimension 2, if a = b = 1, one can write a 
decomposition with A^ 1 /2 terms. 

4. Numerical Results 

We will present in this Section some numerical results for the space homogeneous Boltz- 
mann equation in dimension 2, with Maxwell molecules. We will compare the fast DVM 
method presented in Section [3] with the method introduced in [30J (this latter method 
shall be referred to as the classical DVM one). The time discretization is performed by a 
total variation diminishing second order Runge-Kutta method. 

The first remark concerning the numerical simulations is that, thanks to the discrete 
velocity approach, the conservations of mass, momentum and energy is only affected by 
the aliasing error and thus, for a sufficiently large computational domain, it is exact up 
to machine precision. This is a relevant advantage compared to the spectral (classical of 
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Number of 


Classical 


Fast DVM 


Fast DVM 


Fast DVM 


Fast DVM 


points N 


DVM 


with N = 1 


with N = 3 


with N = 7 


with N = 14 


8 


1.445E-3 


1.4511E-3 


X 


X 


X 


16 


8.912E-4 


9.887E-4 


8.9646E-4 


X 


X 


32 


6.1054E-4 


6.5209E-4 


5.8397E-4 


6.1328E-4 


X 


64 


2.6351E-4 


4.094E-4 


2.906E-4 


3.667E-4 


2.7341E-4 


128 


X 


2.6669E-4 


1.8245E-4 


2.0371E-4 


1.6341E-4 



Table 1 . Comparison of the L error between the classical DVM method 
and the fast DVM method with different values of N at time T = 0.01, 
after one iteration. 



fast) methods, where only mass (and momentum if one considers symmetric distributions) 
is conserved exactly. 

Let us now present some accuracy tests. In the case of two dimensional Maxwell 
molecules, we have an exact solution of the homogeneous Boltzmann equation given by 



f(t,v) 



exp(-v 2 /2S) 



2S-1 + 



1-5 
2S 



2ttS 2 

with S = S(t) = 1 — exp(— 1/8)/2. It corresponds to the well known "BKW" solution, 
obtained independently in [2] and |23| . This test is performed to check the accuracy of 
the method, by comparing the error at a given time T en d when using N = 8 to N = 128 
grid points for each coordinate (the case N = 128 for the classical DVM has been omitted 
due to its large computational cost). We give the results obtained by the classical DVM 
method and the fast one, with different numbers of N. We choose the value N such that 
the classical method is convergent according to Theorem 12.41 namely 

2N 



N 



_3 + \/2j ' 

Then, one has N = 1 when N = 8, N = 3 when N = 16, N = 7 when N = 32 and 
N = 14 when N = 64. These values give a result corresponding to the kernel mode (|3.2|) ; 
namely that no truncation of the number of lines has been done: the solution obtained 
is essentially the same obtained with the classical DVM method. Note that N must be 
chosen less of equal than N and this is why we do not present the results with, e.g., N = 16 
and N = 7. 

Table Q] shows the relative L 1 error between the exact "BKW" solution / and the 
approximate one /j. It is defined by 

El- N \fi(t)-f(vi,t)\ 



S 1 (t) 



£f=-AT 1/^)1 

The size of the domain has to be chosen carefully in order to minimize the aliasing error. 
In this test, we used T = 5 for N = 8, T = 5.5 for N = 16, T = 7 for N = 32 and T = 8 
for N = 64, 128. 
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We can see that, even with very few directions, there is a small loss of accuracy for the 
fast DVM method compared to the classical one, and that taking all possible directions 
we recover the original DVM solution. The observed order of convergence in N is close 
to 1, as predicted by Theorem 12,41 and nearly the same for all values of the truncation 
parameter N (with a small loss for N = 1). 

We also observe that the method is convergent with respect to N, although being not 
necessarily monotone (in the sense that the accuracy can be better for a fixed couple 
of parameters (JV, Ni), Ni < N, compared to the result obtained with another couple 
(iV, N2) with N\ < N2 < N). This is due to the very irregular discrete sphere associated 
with the Farey series, which implies that the information contained in the kernel modes 
can be more complete with the Farey series J 7 ^ rather than • 

We then compare in Figure [5] the time evolution of this error, still in L l norm. We can 
see that it increases initially (exactly as in the classical and fast spectral methods |17|). 
and then decreases monotonically in time. A saturation phenomenon due to aliasing errors 
finally occurs as for the fast spectral method (see [T7], Figure 1). 
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Number of 
points N 


Classical 
DVM 


Fast DVM 
with N = 3 


Fast DVM 
with N = 7 


Fast DVM 
with N = 14 


Fast DVM 
with N = 28 


16 


2 s. 95 


s. 5 


X 


X 


X 


32 


2 min. 18 s. 


3 s. 19 


14 s. 52 


X 


X 


64 


133 min. 44 s. 


16 s. 2 


73 s. 4 


4 min. 43 s. 


X 


128 


X 


85 s. 8 


6 min. 18 s. 


23 min. 2 s. 


92 min. 11 s. 



Table 2. Comparison of the computational time between the classical 
DVM method and the fast DVM method with different values of N at time 
T= 1 {At = 0.01). 



We then give the computational cost of the classical and fast DVM methods in Table 
[2J Here one can see the drastic improvement when comparing the two methods: taking 
e.g. N = 64 points in each direction, the fast method is more than 28 time faster than 
the classical one when no truncation is done (i.e. when we take N = N = 14), and even 
109 times faster with a small loss of accuracy when taking N = 7. 

We also present the evolution of these computational times with respect to the total 
number of points in Figure It is clear when we look at the interpolant curve that 
the theoretical predictions and the effective computational costs agree perfectly. When 
N is fixed, the fast DVM method is of order iV 2 log(iV) whereas when N is fixed, the 
dependence in N is very close to N 2 (actually, the slope of the interpolant curve is about 
1.9). 

5. Conclusions 

We have presented a deterministic way for computing the Boltzmann collision operator 
with fast algorithms. The method is based on a Carleman-like representation of the 
operator that allows to express it as a combination of convolutions (this is trivially true 
for the loss part but it is not trivial for the gain part). A suitable periodized truncation of 
the operator is then used to derive fast algorithms for computing discrete velocity models 
(DVM). This can be adapted to any DVM, provided it features a decoupling properties 
on the quadrature nodes. Our approach will bring the overall cost in dimension d to 
0{N d N d log 2 N) where N is the size of the velocity grid and N is the size of the grid used 
to compute directions in the approximation of the discrete operator. Numerical evidences 
show that the quantity N can be taken small compared to N. Consistency and accuracy 
of the proposed schemes are also presented, both theoretically and numerically. 

Acknowledgments. The first author wishes to thank Bruno Sevennec for fruitful dis- 
cussions on the Farey series. The third author wishes to thank Francis Filbet for fruitful 
discussions and comments about the implementation of the numerical method. 
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Total number of grid points n := A^ 2 

Figure 3. Evolution of the computational time with respect to the total 
number of points for the classical and fast DVM methods, at time T = 1 
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